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Abstract 

We simulate lattice QCD at non-zero baryon density and zero temperature 
in the quenched approximation, both in the scaling region and in the infinite 
coupling limit. We investigate the nature of the forbidden region - the range 
of chemical potential where the simulations grow prohibitively expensive, and 
the results, when available, are puzzling if not unphysical. At weak coupling 
we have explored the sensitivity of these pathologies to the lattice size, and 
found that using a large lattice (64 x 16 3 ) does not remove them. The effective 
potential sheds considerable light on the problems in the simulations, and 
gives a clear interpretation of the forbidden region. The strong coupling 
simulations were particularly illuminating on this point. 
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I. INTRODUCTION 



The absence of a simulation algorithm for finite density QCD is an outstanding problem 
both of QCD thermodynamics and lattice field theory [1]. From the lattice viewpoint, the 
source of these difficulties is the complex nature of the QCD finite density Action [2], [3]. 
This property prevents the use of naive probabilistic methods in evaluating the functional 
integral, thus calling either for exact studies, which are extremely compute intensive, or for 
suitable approximations. 

Many studies, especially the early ones [2], [4], have pursued approximations. Finite 
density QCD has been studied in the quenched approximation, which is the most widely 
exploited approximation in lattice QCD: this approximation has problems in the chiral limit 
because it ignores the axial anomaly, but it is satisfactory for many practical calculations. It 
was therefore rather perplexing and disappointing when severe problems were first reported 
in quenched simulations of finite density QCD [4]: at finite quark mass, the onset of chiral 
symmetry restoration appeared to occur at a chemical potential of half the pion mass, 
which would extrapolate to zero in the chiral limit. Several explanations of this seemingly 
unphysical result were proposed, but none of them led to a satisfactory treatment of finite 
density QCD, nor did they indicate if the quenched approximation was at fault (except in 
that they relied on simple probabilistic arguments that failed for the complex action of full 
QCD at finite chemical potential) or if there were intrinsic problems in the lattice formulation 
of fermions and chemical potential that would survive more sophisticated treatments. 

In Ref. [5], we attempted to confirm the relationship fi onset = in the quenched 

model. We considered that to be necessary because, as suggested by other authors in the 
past, the coincidence of the onset with one half the pion mass might only have been a numer- 
ical accident, the correct relationship being \x onse t = — A, where A is the contribution 
of the nuclear binding energy. If this was the case, the problems with finite density would 
not have been serious. These considerations called for high precision measurements, never 
performed before. Later arguments [6] indicating that the transition should indeed occur 
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precisely at \x = lack rigor and could well do with the scrutiny of such measurements. 

This work, alas, confirmed the onset at fi = m v /2. 

However, as discussed at length in Ref. [5], the finite temporal extent of the lattice 
might be responsible for unphysical closed quark loops winding around the system. Such 
loops can generate huge fluctuations and give misleading indications of deconfinement and 
chiral symmetry restoration. These considerations motivated us to follow-up Ref. [5] with 
simulations on a lattice of greater temporal extent. We discuss these results in Section II 
below. In particular, we shall give an overview of the results in Section II. A, while Sections 

II. B and II. C , are devoted to the spectrum analysis and finite size effects, and can be 
skipped by the reader uninterested in technical details. We next turn to the strong coupling 
limit of the theory. Section III is devoted to new simulations performed in that regime. The 
strong coupling limit has little phenomenological relevance, but it displays both confinement 
and chiral symmetry breaking - the QCD characteristics which are relevant in this context. 
There are several considerations which made this study appealing to us: first, in the same 
spirit of the weak coupling simulations, we wanted to have a clean measure of the onset: at 
strong coupling half the pion mass and one third the baryon mass are far apart, and are 
easily distinguished numerically; second, the simulations are much less computer intensive, 
and accurate results can be obtained; finally, several analytic treatments are available, which 
can clarify the interpretation of the results. We have also introduced new observables which 
have helped clarify some pathologies, and have helped us arrive at a coherent description 
of the failure of the quenched approximation. We discuss these measurements in Section 

III. A. Section III.B contains the ordinary thermodynamics results. They were particularly 
clear, and were easily extended to [i > niN/3, the saturation region. Finally, in Section III.C 
we contrast the results of the simulations with the analytic predictions. This gives a new 
interpretation of the forbidden region, \x > 
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II. QCD AT NON-ZERO BARYON DENSITY IN THE SCALING REGION: 
IMPROVED OPERATORS, TEMPERATURE EFFECTS AND WINDING 

LOOPS. 

This section contains an upgrade of our previous work on quenched QCD at non-zero 
baryon-number density and j3 = 6.0 [5]. We briefly summarize here our main findings and 
open problems, and refer to [5] for introductory material and a general review. Our best 
success there was a measurement of the nucleon mass for \x < m n /2. The results were in 
agreement with the predictions of a simple constituent quark model, and hinted at a critical 
value for the chemical potential given by one third the baryon mass. Unfortunately, the more 
interesting region fi > evaded us: enormous fluctuations dominated the results, and 

not only was it impossible to estimate the mass spectrum, but the shape of the propagators 
themselves were grossly distorted. The results for thermodynamic quantities were somewhat 
better, and we were able to detect a modest decrease in the chiral condensate at fi ~ m n /2. 

The main purpose of this new work is thus twofold: first, to understand and possibly 
eliminate the source of the anomalous fluctuations in the spectrum, in order to extend our 
simulations inside the interesting region of the chemical potential; second, to explore the 
possibility of residual temperature (the temperature being the reciprocal of the time extent) 
effects on the thermodynamic observables. Both of the above called for a simulation on 
a lattice of larger temporal extent: this pushes the temperature further toward zero, and 
possibly suppresses precocious, unphysical quark loops winding around the lattice in the 
time direction. We have also proposed another possible explanation for the pathologies 
afflicting the spectrum: the wall source used for the spectrum measurements could produce 
strong distortions in the propagators. Remember that such a source, used in conjunction 
with with a point sink, alters the value of the coefficients aj in the expansion of a propagator 
G(t), G(t) = J2 a i emit , by maximizing ao, the projection onto the fundamental state, a clear 
numerical advantage. Unfortunately, this advantage is accompanied by a possible lack of 
positivity of the other coefficients, leading to instabilities in the results, amplified by the 
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asymmetries introduced by the chemical potential. In the new set of simulations we have 
thus decided to use a noisy-wall source: namely, a random SU(3) matrix was put on the 
(even, even, even) sites of the spatial cube t = , and the Dirac operator was inverted in 
this background. It is easy to show that the average local-hadron propagators are the same 
as the ones obtained with a point source, thus avoiding the positivity problems we have 
just mentioned. There is a statistical gain due to the increase of the number of source 
points (roughly, its square root), and a loss because of possible incomplete cancelations of 
gauge- variant terms which are inherent in this technique [9]. The real advantages can only 
be judged a posteriori. A sure disadvantage, well known from the zero chemical potential 
simulations, is the absence of useful results for the baryon - with this source, asymptotic 
behavior sets in only at a large time interval where the baryon propagator is too small to 
measure. (N.B. we considered using a version of the noisy-wall source which corresponds 
to smeared sources, along with the corresponding smeared sinks. However, this gives less 
statistical gain than for the point version.) 

Summarizing, we have doubled the temporal size of the lattice ( 16 3 x 64 vs. 16 3 x 32), 
and used a noisy-wall (as opposed to a rigid wall) source for the spectrum measurements. 
Making these two changes (source, and time extension) simultaneously was the best way 
to optimize our computer resources, but unfortunately we were not able to disentangle the 
main sources of the troubles afflicting previous spectrum computations. The measurement 
strategy for the order parameter and number density was the same as before, so in these 
cases we are genuinely monitoring finite temperature effects. 

As in past work, we exploited the "global" Z 3 symmetry of quenched QCD, using three 
different anti-periodic boundary conditions for the fermions in the time direction: tpi(N t ) = 
—Ziip(0), where Zi stands for each of the three cube roots of unity. Remember that this can 
help enforce the constraints of confinement by suppressing unphysical winding loops. It is 
thus interesting to monitor the statistical gain coming from averaging over these boundary 
conditions as a function of the lattice time extent. The bigger the gain, the more important 
the role of the (unphysical) winding loops, so the interplay of temperature and winding loops 
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can be clearly exposed. 

We have generated 30 independent configurations on a 16 3 x 64 lattice from a code which 
is a blend of Metropolis and over-relaxed algorithms. A subset of these configurations was 
used for the inversion of the Dirac operator for each n value. The quark mass was .02, as 
in our previous work. We have inverted the Dirac operator for fi = 0.0, .1, .15, .17, .2 on 
13, 13, 26, 12, 28 configurations, respectively. We have also done a few exploratory runs for 
larger fi values (fi ~ 1.), searching for hints of saturation (see [4], and the following Section), 
but the number of iterations required for the inversion was always huge, and the number 
density rather low. In the free theory saturation occurs only for rather large \x : we probably 
observed a similar trend at j3 = 6.0, in agreement with the observations of [16]. 

The measurements and analyses were exactly the same as in our past work, to which we 
refer the reader for details. 

A. Overview of the results 

The results for the chiral condensate and the number density are shown in Table I and 
II and collected in Figs. 1 and 2. They are systematically compared with those obtained 
on the smaller lattice (for the sake of comparison, we have repeated the analysis on the 
smaller lattice including all the data, without any cuts, which accounts for some difference 
(statistically irrelevant) among the results in this table and the ones of Ref. [17]). 

For small values of the chemical potential (// < we do not see any temperature 

effect for the chiral condensate, while for // > m v /2 the overall behavior does not exclude 
the possibility of (modest) finite temperature effects. We will come back to this point in the 
last subsection. 

In the spectrum the most visible phenomenon concerns the pion propagator. The huge 
and unphysical fluctuations described above in our past measurements [17] disappear. How- 
ever, the number of iterations required for the inversion when [i exceeds m v /2 remains big, 
of order 10 4 , meaning that the Dirac operator is still nearly singular. 
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< ipip > 




16 3 x 32 


16 3 x 64 


0.0 


.1377(6) 


.1376(2) 


.1 


.1375(8) 


.1376(3) 


.15 


.1362(18) 


.1381(5) 


.17 


1359(1401 


1357(191 


.20 


.1088(100) 


.1277(35) 




TABLE I. Results for the chiral condensate as a function 


of the chemical potential on the 


16 3 


x 32 and on the 16 3 x 64 lattices. We quote the average over 


the three boundaries and the two 



opposite jj, values 



< Jo > 

H 16 3 x 32 16 3 x 64 

0.0 .00016(93) .00029(50) 

.1 .00033(83) .00011(56) 

.15 .00069(145) -.00013(52) 

.17 .00071(161) .00028(212) 

.20 .00510(899) .00481(315) 



TABLE II. Results for the number density as a function of the chemical potential on the 16 3 x 32 
and on the 16 3 x 64 lattices. We quote the average over the three boundaries and, for \i ^ 0, over 
the two opposite \i values 
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Figure 1. Results for the chiral condensate at (5 = 6.0 (diamonds). The results on the 16 3 x32 
lattice are shown for comparison (squares) 
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Figure 2. Results for the number density at (5 = 6.0 (diamonds). The results on the 16 3 x 32 
are shown as squares, slightly displaced. 
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If we ascribe the improvement to the different sources we were using, we would conclude 
that the source affects the amplitudes, but not the poles, as expected. Despite the better 
behavior of the propagators, the pion mass estimate was not compelling for \x > m n /2, while 
for smaller chemical potentials we confirm previous results. We shall give some detail on 
the mass analysis in the next subsection. Here, we quote the results for the pion mass in 
Table III, again comparing with the results on the smaller lattice. We see that there are no 
discernible temperature effects. Also, remember that we were using two different sources, 
wall, and point with the noisy improvement. For our wall source, the effective mass estimator 
appears to reach its asymptotic behavior from below (see the details in the analysis Section 
below). In particular, as we will show in more detail later, our mass estimate at [i = .17 is 
certainly non-asymptotic and should be considered as a lower bound. On the other hand, 
because of the positivity properties we mentioned above, the effective mass estimates from 
the propagators obtained with noisy-wall source give upper bounds. One might be tempted 
to "average" the results from the two different sources, thus obtaining a nearly constant 
pion mass across the transition. 

We summarize all our findings for the spectrum in Fig. 3. Besides the old results, and 
the new ones from Table III, we show there the estimates for the pion mass at \i = .2 (see 
next section for discussions and caveats), the point with the huge error coming from the 
full sample, the other one, slightly displaced, from a subset from which the most wildly 
fluctuating propagators have been excluded (in the preliminary results reported in [17] the 
statistical sample was smaller, so the cut was mandatory and we did not have an estimate 
from the full configuration ensemble). The fairest statement is that, despite our attempts, 
we do not have a reasonable mass estimate at fi = .2, even though the propagator was better 
behaved than the one on the smaller lattice. 

It would be interesting to determine if the pion mass does indeed increase for \i > 
which would be an indication of chiral symmetry restoration. Alternatively, we can search 
for hints of chiral symmetry restoration in the spectrum sector by looking at the propagators 
of the chiral partners. We have indeed observed (Fig. 4) that the scalar and pseudoscalar 
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H 16 3 x 32 16 3 x 64 

0.0 .3396(36) .3379(31) 

.1 .3374(75) .3440(51) 

.15 .3182(52) .3464(142) 

.17 .313(15) -41(3)(4) 



TABLE III. Pion masses as a function of \i on the two lattices. See text for details 





0.0 


.1 


.15 


.17 


.20 


.20 (cut) 


F2 [5-32] 


.3379(31) 


.3440(51) 


.3464(142) 








Em [20-30] 


.3356(42) 


.3376(47) 










F2 [2-9] 


.3610(55) 


.3510(55) 


.3609(57) 


.4107(324) 


.58(29) 


.32(13) 


F2 [3-7] 


.3913(88) 


.3914(87) 


.3807(104) 


.4155(433) 


.60(24) 


.37(24) 


Em [14-16] 


.3440(45) 


.3442(52) 


.3307(584) 









TABLE IV. Estimates for the pion mass from effective masses (Em) or two particle fits (F2). 
The first two lines are the "good" results 
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gure 3. Summary of the spectroscopic results on the two lattices. See text for details. 
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Figure 4. Scalar (circles) and pseudoscalar (crosses) propagators as a function of \i. The 
(upper, middle, lower) couple of values is for the propagators at time separation (0,1,2). 
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propagators tend toward degeneracy. The same holds in the vector-pseudovector channel, 
but in both cases the effect is not dramatic, given that the two propagators are very close 
even at \i = 0. Moreover, because of the different sources we were using on the smaller 
lattice, a finite size study of these quantities is not possible at this time. 

In conclusion, we have observed only very modest hints of chiral symmetry restoration in 
the spectrum and in the thermodynamics, possibly affected by residual temperature effects. 
We have shown how suitable observables can help control major pathologies. However, these 
results alone do not determine the physical nature of the onset at [i = m n /2. 

B. Details of the spectrum analysis 

We present here some details of our spectrum analysis. They provide some information 
on the qualitative effects of the chemical potential on the propagators. 

We concentrate on the pion propagator, which is the best-behaved. The baryon propaga- 
tor loses its numerical significance at distance 1-2 and we cannot extract useful information 
from it. 

We show in Fig. 5 the collection of the effective mass plots for the pion (later in this work 
we will also extract masses from the t-dependence of the propagators). The quality of the 
data is good for small fi, and plots, such as Fig. 6 for // = .1, show good agreement between 
our mass estimates obtained from the two lattice sizes, and the two source methods. As 
[i is increased, the signal is lost earlier and earlier in time. A better idea of the quality of 
the results can be obtained from Table IV, where we give the mass estimates from different 
analyses. It is clear that the asymptotic behavior for the effective mass sets is only for t > 20 
(considering the behavior of fi = .0 and fi = .1), so the effective mass analysis does not give 
results for /j, > .1. 

Next we considered two particle fits (fundamental + excited). In this case we were able 
to use time separations as small as five, although large time distances were still needed to 
stabilize the results. Fig. 7 demonstrates the quality of the fits, which were successful up to 
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Figure 5. Effective masses as a function of time for n = (0.0, .1, .15, .17, .2) 
(pluses, circles, squares, diamonds, crosses). Data for different /i's are slightly displaced. 
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Figure 6. The results for the effective mass at /j, = .1 on the big lattice are contrasted with 
the ones on the smaller lattices, for the three different boundaries. 
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Figure 7. Two particle fits to the pion propagator in the time interval [5-32], at n = .1 
(circles) and .15 (crosses). The fitting curves are practically coincident. Only the data 
points statistically different from zero are shown. 
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Figure 8. Two particle fit to the pion propagator at \i = .2 in the interval [0-9] (solid line). 
A similar fit, but with the mass in the fundamental channel constrained to its fi = value 
is shown as a dashed line. 
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To get mass estimates at /j, = .17 and \x = .2 we first restricted the fitted interval to a 
statistically significant range (1 < t < 10). Then we estimated the error due to the lack 
of asymptotic behavior by comparing, at \i = 0, the result obtained on that range with 
the correct one, and we assumed that such an estimate is \x independent. In this way we 
estimated the mass and errors at fi = .17, given in Table III, the first error being statistical 
and the second computed as just said. It is clear that this mass is certainly not asymptotic, 
and has to be considered only 'indicative'. The same problems occur at fi = .2, where we 
have also tried some fits which discard the noisiest propagators. The situation at fi = .2 is 
best described by a figure. We show in Fig. 8 two fits to the pion propagator at \i = .2. The 
solid line is drawn in correspondence to the unconstrained fit, which gives for the central 
mass m n = .60, with a statistical error A(m 7r ) = .24. The dashed line is for a fit with a 
constrained value in the fundamental channel, m v = .3379 as in the \x = 0.0 case. The 
two fits are coincident in the region where the propagator is statistically different from zero. 
Only better results at large t would allow a safe estimate for the pion mass. 

C. Further search for finite temperature effects 

We have said that the results for the chiral condensate and the number density do not 
exclude the possibility of finite temperature effects. We have then selected three [i values, 
.0, .15, .20, where we had the same number of gauge configurations on the two lattices in 
order to do a more detailed comparison. 

In Table V we record the results for the chiral condensate for the different boundary 
conditions, and opposite \i signs, in order to find systematic trends in their average values, 
or in their fluctuations, [i = sets the scale for the expected effects. We do not see any 
size effects at [x = 0, while the error ratio is roughly V2, expected on a purely statistical 
basis (the lattice volumes differ by a factor of 2). Deviations from this trend have to be 
interpreted as induced by the chemical potential. At [i = .20 it appears that the difference 
in the results is systematic, and statistically significant. 
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.1145(90) 


.1230(154) 


.0843(248) 


.1301(48) 


.1328(82) 
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.1292(55) 


.1073(142) 


.0944(267) 


.1262(72) 


.1250(43) 


.1237(36) 


TABLE V. Temperature dependence of < tptp > . 


We compare the results for the 16 3 x 32 and 


16 3 x 64 lattices at \x 


= 0, +.15, +.20, for the three different boundary conditions. 




<Jo> 






16 3 x 32 






16 3 x 64 








z 2 


z 3 


Zi 


z 2 


^3 


0.0 


-.00048(99) 


.00081(117) 


.00063(128) 


.00034(88) 


.00039(88) 


.00015(90) 


+.15 


-.00040(123) 


.00228(297) 


-.00219(140) 


.00024(91) 


.00170(72) 


.00081(101) 


-.15 


-.00232(123) 


.00291(297) 


.00078(161) 


.00194(98) 


.00053(83) 


.00107(95) 


+.20 


-.01665(973) 


-.00311(1350) 


.04192(3721) 


.00548(499) 


.00901(585) 


.00819(541) 


-.20 


-.01253(736) 


-.00023(731) 


.00433(1263) 


-.00886(447) 


-.00207(777) 


.00477(324) 



TABLE VI. Temperature dependence of Jo . We compare the results for the two lattices at fx 
= 0, +.15, +.20, for the three different boundary conditions. 
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Figure 9. Collection of the results for the number density at (5 = 6.0, for the three different 
boundaries, as a function of //. The results for // and —fi have been averaged. The circles 
are for the small lattice, the squares, displaced, for the big one. 
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Also, the error ratio grows larger, thus indicating some non-trivial improvement mechanism 
on the larger lattice, probably induced by the suppression of winding loops. Note also that 
the statistical gain coming from averaging over the boundary conditions is bigger on the 
smaller lattice, and that on the smaller lattice there is a much bigger dispersion from one 
boundary condition to another. This suggests that winding loops are playing a significant 
role, which diminishes as the lattice size increases. 

Table VI, where the results for the number density are given, shows that J is always 
consistent with 0. From Fig. 9 we see that on the smaller lattice (circles) there are exceptional 
events (deconfined configurations) which are absent on the bigger lattice (squares). Note 
that some of the exceptional configurations have a density of the "wrong" sign: for instance, 
a negative chemical potential is apparently inducing a positive density! 

In conclusion, by increasing the temporal extent, lattice artifacts are lessened. However, 
the improvements are very modest, so there is no reason to believe that the problems inherent 
in quenched simulations will be solved simply by using large volumes. 

III. THE INFINITE COUPLING LIMIT OF NON-ZERO DENSITY QCD 

Recall, once more, the nature of the finite-density quenched QCD puzzle : observables 
deviate from their // = 0.0 values when // = m v /2. In particular, the chiral condensate 
appears to fall at \i = m n /2. This result, which has suggested to other workers in the field 
that an infinitesimal value of the chemical potential restores chiral symmetry in the chiral 
limit, has been observed in simulations in the scaling region, as well as in the strong coupling 
limit [4]. At strong coupling this onset of sensitivity to \i observed in numerical simulations 
also differs from analytic predictions. This is rather peculiar, since analytic mean-field 
strong-coupling results usually compare well with numerical simulations. It will prove to 
be informative to uncover the reasons behind the discrepancy at strong coupling because 
it will suggest the reason for the failure of the quenched approximation in general. These 
considerations, together with the low computational costs of these simulations (compared 
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to the simulations in the scaling region described above), motivated the following study. 

We generated 20 random configurations on a 8 3 x 16 lattice, and 10 configurations on 
a 8 3 x 32 lattice. Again, we used three different anti-periodic boundary conditions for the 
fermions in the time directions : ipi(t) = — Zj-0(O), where Zi stands for each of the three 
cube roots of unity. We first chose a bare quark mass of .1. The strong coupling predictions 
for the pion and baryon masses at m q = .1 are .6 and 3.3, respectively. In our simulations \i 
ranged from 0. to 1.2, thus including the interesting region .3 < fi < 1.1 and the mean field 
prediction for the pseudo-critical [i ( [i c ~ 0.6). Moreover, a subset of 7 configurations on the 
smaller lattice was analyzed with a very heavy quark mass (m q = 1.5), with \x in the interval 
(1.2:1.5), the pseudocritical point in this case being [i c ~ 1.37. Finally, 10 configurations on 
an 8 3 x 32 lattice with m q = .1 were used to monitor possible temperature effects. 

We have measured the chiral condensate, the number density, the energy density and the 
pion mass. The other masses are too heavy, and difficult to extract - a typical limitation 
of strong coupling simulations. Along with standard observables, the unphysical operator 
G fl G^ (the "false" or "baryonic" pion, suggested by [6] as opposed to the real pion G^Gt^) 
was measured to obtain information about the poles of the Dirac operator. 

The main results of this study are : 

First, there are three distinct intervals of chemical potential. A conventional region, for 
\i < m n /2, where all the zero-triality observables maintain their \i = value. A "forbidden" 
region, for m v /2 < [i < m^/3 (note the dependence on the nucleon mass), characterized by 
wild fluctuations, and by average values which (when statistically meaningful) show devia- 
tions from their // = 0.0 values. Finally, a saturation region, where all the thermodynamic 
observables have their limiting, large [i values, the amplitudes in the mass spectrum drop 
to zero, and no single particle states exist. 

Moreover, inspection of the effective potential as computed in a mean field approximation 
shows that the "forbidden" region, which begins at fi = corresponds to the onset of 

metastability in the effective potential. 
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Figure 10 . The propagator of the "baryonic" pion is shown as a function of time for /j, = 
0.0, .1, .2, .3, .4, .5, .9. Note the T-asymmetry for fi < m v /2 , and the flatteness associated 
with the forbidden region. 
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Figure 11. The propagator of the "baryonic" pion for \x = 1.0, 1.1, 1.2, in the saturation 
region. 
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gure 12. "Pion" effective mass versus time, at \i = (0.0, .1, .2, .3), from top to bottom 
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Figure 13. "Pion" masses, from an effective mass analysis, versus the chemical potential \x. 
The dashed line is y = m^O) - 2/j. 
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A. Fluctuations, poles of the Dirac operator and the "baryonic" pion 

As discussed at length in past work, when the chemical potential fi exceeds half the pion 
mass, the fluctuations increase dramatically. This can be read off the eigenvalue distribution 
[2], since when n ~ m v /2 the quark mass lies inside the eigenvalue cloud, thus producing 
eigenvalues of the full Dirac operator close to zero. These huge fluctuations are associated 
with a dramatic increase in the CPU time required for the inversion of the Dirac operator: 
the average number of conjugate gradient iterations for \x < m^/2 was ~ 100, while when \i 
exceeds m n /2 the number of iterations grows large (~ 10 4 ). Again, this is consistent with 
the eigenvalue picture, the number of iterations being proportional to the ratio of the largest 
to the smallest eigenvalue. For fi > 1 the number of iterations is again very small (~ 10). In 
this region we are inside the eigenvalue "crescent" , and the first Brillouin zone is saturated. 
Again, we note that this value corresponds, approximatively, to one third the baryon mass: 
the theory, at least in its strong coupling limit, is sensitive to the nucleon, as predicted by 
simple nuclear models. 

In order to gain additional insight into the eigenvalue distribution and zero modes of the 
Dirac operator, we studied a suitable, "unphysical" operator, the "baryonic pion" (which 
we shall often denote "7r"): G^G^, as opposed to the real pion G^Gt^, G^ being the quark 
propagator computed with chemical potential //. We shall see that the amplitude of the 
baryonic pion propagator grows very large when \x > m^/2, and its poles, which are related 
to the poles of the Dirac operator, are very well exposed. 

We show the 'V propagator in Figs. 10 and 11. For small values of the chemical potential 
the propagator has a baryon-like (T- asymmetric) behavior . For \x > m^/2 the propagator 
flattens, which can only happen because of zero eigenvalues of the Dirac operator: this gives 
direct evidence of the zero modes we were looking for (Fig. 10). In the saturation region 
(Fig. 11), G" 7r "(t) is again baryon-like, but this time the propagation, as expected, is in 
the opposite direction. Outside the forbidden region, which, in this case, is associated with 
the flat propagators, we can measure the baryonic pion mass. For that purpose, we relied 
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on a simple effective mass analysis, and evaluated the logarithmic derivative without any 
attempt at symmetrization (we hesitated to use the same parameterization as for staggered 
baryons). The flat region in the effective mass plots extends till the influence of the backward 
propagating state becomes appreciable : from our results it is clear that this still produces 
a reasonable interval, so that we can obtain a safe estimate for the baryonic pion mass. 
In Fig. 12 we show the effective mass plots, together with the results of a conventional 
effective mass analysis for the real pion. The 'V mass satisfies, for // < m v /2, m<v' (/-*) = 
m<v' (0) — 2//, (where m<v> (0) = m 7r (0) = m v ) thus extrapolating to zero for // = m^/2 (Fig. 
13). 

In conclusion, the study of the baryonic pion gives clear evidence of zero modes in the 
quark propagator for \x ~ m^/2, in agreement with the eigenvalue picture. 

Summarizing, these measurements provide a coherent description of the spectral struc- 
ture of the Dirac operator over the entire fi range. The poles for fi > m v /2 can be clearly ex- 
posed by suitable, yet unphysical observables. Their mathematical and physical significance 
is open to many interpretations which have been discussed and reviewed in the literature. 
We shall comment further on this point below. 

B. Physical observables 

We have measured the particle spectrum, and thermodynamic observables (chiral con- 
densate and number density), via a standard stochastic estimator. 

The pion propagator (Fig. 14) is completely insensitive to the chemical potential up to 
m n /2. The "forbidden" region (m v /2 < \i < m B /3) is dominated by fluctuations, while in 
the saturation region (/i > m^/3) the pion propagator amplitude drops to zero. We were not 
able to measure the baryon mass (a common drawback of strong coupling simulations), thus 
missing a very important piece of information about the critical behavior. It was however 
possible to monitor the amplitude of the baryon propagator. At m ? = .1 we found that it 
decreased by an order of magnitude in the saturation region ( the amplitude is ~ .01 in the 
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Figure 14. Pion propagators at n = (0.0, .1, .2, .3) (Plus's, crosses, diamonds, circles.) The 
pion mass is insensitive to the chemical potential in this range of values. 
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conventional region, undefined in the forbidden one, and lO^" 6 ' 9 ' 11 ) at = (1, 1.1, 1.2)). This 
trend in the amplitude is common to all the propagators we have measured, and indicates 
that no real particle can exist in the saturation region. 

The results for the thermodynamic observables are shown in Figs. 15, for m q = .1 
(Fig. 15a) and m q = 1.5 (Fig. 15b). As noticed in past work, the chiral condensate and 
the number density are closely correlated. Saturation is especially clear in the number 
density, which is three -the maximum number of quarks which can occupy a single site- 
for \x > mjv/3. Our results are consistent with those previously reported in the literature, 
whenever a comparison is possible (fi < m n /2). 

It is also informative to consider the results configuration by configuration. We know 
that ensemble averages have some special features at finite density, and a more detailed 
analysis can be worthwhile. For that purpose, we use the data at m, = .1 which has the 
best statistics. 

In Figs. 16 we show the number density as a function of the configuration number. At low 
fi (Fig. 16a), the densities computed with opposite values of the chemical potential (joined 
with dotted and dashed lines, respectively) are strongly correlated: so we do not see any 
sensitivity to the chemical potential itself. In the intermediate region (Fig. 16b), J is wildly 
oscillating. We observe configurations which have a non-zero density, but with the wrong 
sign, for instance when a positive density is obtained with a negative fi. This, once more, 
tells us that one must be careful in using concepts like confinement and chiral symmetry 
breaking on an isolated configuration: a negative chemical potential should enhance, in the 
statistical average, antibaryon propagation, while we see that on isolated configurations the 
opposite can happen. In the saturated region (Fig. 16c) things are very clear, and each 
configuration has a net density with the "right" sign. 

In Figs. 17 and 18 we give the full ensemble of values for the chiral condensate and the 
number density, as a function of the chemical potential. For the number density, we took the 
average over ±/i, which enforces T-symmetry event by event, and produces J = at fi = 0. 
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Figure 15a. Summary of the thermodynamics results at strong coupling. The chiral con- 
densate ( circles, left ) and the number density ( crosses, right ) are plotted as a function of 
[i , for m q = .1 
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Figure 15b. As in Fig. 15a, but m q = 1.5 
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Figure 16a. The results for the number density JO are shown configuration by configuration, 
for the three boundaries and the two opposite fi values at [i = .1 
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Figure 16b. As in Fig. 16a, but \x = .5 
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Figure 16c. As in Fig. 16a, but \x = 1.2 
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Figure 17. Collection of the results for the chiral condensate, for the three differ- 
ent boundaries, and the opposite values of function of the chemical potential. 
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Figure 18. Collection of the results for the number density, for the three different boundaries, 
as a function of //. The results for fi and —ji have been averaged 
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For n > m n /2 the results for < tptp > and J obtained configuration by configuration 
spread over a wide range, which includes both the saturation and the [i = values. 

The behavior we have described suggests that we might have entered a metastable region 
rather than the symmetric phase, and that the identification of the onset of \x dependence 
in the theory's observables with the pseudocritical point computed in a mean field approxi- 
mation may not be correct. In the next subsection we shall compare this interpretation of 
the data with the predictions of the strong coupling analytic calculations. 

C. Effective models and simulation results 

Several, qualitatively equivalent, strong coupling, 1/d expansions of full QCD with stag- 
gered fermions at finite density and temperature are available in the literature [4] [7] [8] [10]. 
Such mean field analyses predict a strong first order transition at finite chemical potential 
and zero temperature. Since fermion loops do not play a significant role in the leading 
order of the 1/d expansion at strong coupling, it is sensible to try to interpret quenched 
simulations in terms of the analytic mean field predictions. This has been done in the past 
[4] [7], and we will continue in this tradition in some detail. We shall show that the onset 
Ji of the metastable region of the effective potential of these studies is associated with the 
onset fi on set for fi dependent thermodynamics introduced at the beginning of this paper. 

Let us now describe these calculations in more detail. At small \i the free energy, which 
plays the role of an effective potential, plotted as a function of tptp has only one minimum, 
corresponding to the fi = 0.0 value of the chiral condensate < tptp > . While increasing fi 
= p,, the free energy develops a secondary minimum at tptp = . The two minima become 
of equal height for // = //*. // is thus interpreted as the pseudocritical point at finite mass, 
or the exact critical point \x c for chiral symmetry restoration in the chiral limit, {jx — //) 
approaches as the quark mass increases, while these two values of the chemical potential 
are well separated for smaller quark masses, thus creating a wide metastable region. By 
further increasing the chemical potential, the minimum at finite tptp eventually disappears 
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at fi = fi s , the only surviving minimum being the one at ipip ~ 0. This should correspond to 
the complete saturation of the first Briouillin zone, i.e. density = 3. Before proceeding, we 
recall that fi c 7^ m B /2>. This can be understood because of a considerable nuclear binding 
energy at strong coupling [7] (infact, fi c approaches m s /3 when g — > 0), or because of the 
effects suggested in Ref. [16]. Anyway, from the perspective of the present study, the main 
feature of the strong coupling analysis is that the critical chemical potential fi c is different 
from zero, and the lack of exact coincidence with the naive predictions m^/3 carries little 
weight: if the results of the strong coupling analysis were realistic, the theory would not be 
pathological. This is why it is important to reconcile the numerical and analytic results at 
finite mass, or to uncover the reasons behind their differences. 

Summarizing, there are a few interesting values of the chemical potential, labeling the 
different physical regions associated with the chiral transition which emerge from the effective 
potential analysis: fi c , the critical point of the chiral transition in the chiral limit; //, the 
pseudocritical point as analytically computed at finite mass; fi, the onset of the metastable 
region found in an effective potential approach; fi s the chemical potential at which saturation 
occurs. We would like to study the interrelation of these points with fi onset , the onset of 
\i dependence in thermodynamic observables at finite quark mass observed in numerical 
simulations. Naive arguments would suggest that fi* = fi onS et , i-e. that the transition has a 
strong first order character also at finite mass. However, it is easy to imagine situations (as 
we are will explain below) in which fi* and fi onse t are different. In this case, lim m<; ^.o A** = t^c, 
while lJL onse t in itself does not have any immediate relationship with the critical point fi c . 
Rather, it should be identified with Ji : the onset of fi dependence in the thermodynamic 
observables is produced by the edge of the metastability region. We also mention that 
the work of Ref. [13] supports the possibility that fi nset an d fi c are indeed different. In 
Ref. [13] a representation of the partition function at strong coupling written in term of a 
monomer-dimer expansion is used to compute the order parameter and the number density. 
< ipip > deviates from its fi = value at the value fi* of the chemical potential computed in 
a mean field approximation. However the onset fi onse t for fi dependence in thermodynamic 
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quantities, as inferred from the behavior of the number density, occurs much earlier. 

We now turn to a discussion of the results: we shall compare the results for the effective 
potential drawn in Figs. 19a and 19b with the thermodynamic results presented in Figs. 
15a and 15b. 

In the following, we use results and notation from Ref. [7]. We show in Fig. 19a the 
effective potential (after formulae(2.10 and segg. of Ref. [7] ) for the same \i values of our 
simulations at m q = .1, the same quark mass, the asymmetry parameter a T /a s = r = 1, and 
T = 16. Also, we show as a dotted line the effective potential at fi = .7 ~ //*. The curves are 
basically coincident for "large" tptp, while for small values they are qualitatively different. 
In particular, we clearly see the onset of metastabilities at p, ~ .4, and the disappearance of 
the minimum at finite tptp for \x ~ 1. The pseudocritical point // is obscured by the huge 
fluctuations associated with such metastabilities. It is only when /j, > m B /3 ~ /j, s that the 
minimum at finite tptp disappears, and for // > all the observables have their limiting 

values: the chiral condensate is zero and the number density is three. By contrasting this 
behavior with Fig. 15a we see indeed that fi ~ ft, and that saturation occurs for fi > fi s . 
In Fig. 19b we have repeated this plotting exercise for our higher mass value (m q = 1.5), 
with an analogous result. In fact, we have observed that in this case the metastable region 
(Ho < (i < fi s ) in the effective potential shrinks to the expected small interval, as predicted 
from Fig. 15b. 

We learn that the metastable region in the effective potential can be identified with the 
forbidden region observed in the simulations. This clarifies the nature of the phenomena 
observed at m n /2 : the pathologies are not manifestations of an early chiral transition, but 
are indicative of the existence of a region of metastability. For // ~ m v /2 we have the first 
appearance of zero modes. These manifest themselves in huge fluctuations, in the increase of 
the iterations needed for the inversion, and in the flatness of the baryonic pion propagator. 
In particular, the fluctuations observed in the chiral condensate should be associated with 
fluctuations and distortions of the eigenvalue distribution: when the chiral condensate is 
almost zero, the eigenvalue distribution must be such that the point (m q ,0) gets inside. 
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Fig. 19a 



Figure 19a. Free energy, from Ref.7, as a function of ipip, for m q = .1 and the chemi- 
cal potentials used in our simulations, /i increases from top to bottom. The dashed line 
is drawn in approximate correspondence to the mean field prediction for the pseudocritical //. 
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gure 19b. As in Fig. 19a, but m q = 1.5. 
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This, in turn, generates the secondary minimum around zero in the chiral condensate dis- 
tribution, and produces the observed decline in the chiral condensate itself, which thus 
happens much before than the actual pseudocritical point /j,* 1 To master the pathologies for 
\i > m v /2, and expose the physics of the chiral transition which is hidden deep inside the 
metastable region, we need to "enforce" the saddle point solution. 

We thought that a natural possibility would be to increase the time extent of the lattice 
so the system could find the actual minimum : this is why we produced the 10 configurations 
on a 8 3 x 32 lattice mentioned at the beginning of this Section. However, the inversion time 
inside the metastable region remained very large, and the comparison between the results 
of the two lattices excluded strong volume effects. 

This subject deserves further study. Comparisons with the behavior of other models 
can offer some guidance. In particular, we note that the zero temperature, finite density 
transition of the U(l) xU(l) Gross Neveu model, which has been successfully simulated [11], 
is second order [12] , i.e. does not have a metastable region. In SU(N) gauge models the 
transition turns first order for N > 2, while the SU(2) model has a second order transition 
[8]. First order transitions, with their "forbidden" metastable regions, and their numerical 
failures, seem associated with complex actions 2 , while models which can be successfully 
simulated have real actions, and undergo continuous transitions. 

Clearly we need to uncover the physical reasons for the complex action at nonzero yu, 
and/or to explore possible alternatives, such as the Hamiltonian formulation: we note that 



It is possible that the situation improves in the full model, since in that case the eigenvalue dis- 
tribution does depend on the mass value, which can prevent the distribution itself from fluctuating 
randomly around the mass point on the real axes. 

2 This is consistent with the picture sketched above: only an eigenvalue distribution spread over 
the real axes (i.e. generated by a complex action) can produce the secondary minimum at tptp = 0, 
since this is due to the mass point falling inside the eigenvalue distribution. 
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the authors of ref. [18] found that the critical chemical potential, computed at strong coupling 
in Hamiltonian lattice QCD, is indeed equal to the dynamical fermion mass, which, in the 
same scheme, is one third the nucleon mass [19]. 

IV. CONCLUSIONS 

We have performed an exhaustive study of quenched QCD at finite density in the scaling 
region, and in the strong coupling limit. We have measured the standard thermodynamic 
observables, the spectrum, and unphysical, ad-hoc observables meant to elucidate some of 
the peculiarities of the model. We confirmed the pathological behavior observed in the past 
for (i > m v /2. The onset at m v /2 has been measured with great accuracy both at strong 
and weak coupling, and found to be independent of the lattice size, and the measurement 
technique. We think that our simulations in the scaling region have clearly exposed the 
pathologies at m v /2. In particular, it has been shown that the early onset is indeed half 
the pion mass, which rules out the possibility of rescuing quenched, finite density QCD 
by using some refined nuclear matter models. However the accurate simulations performed 
in the scaling region do not suffice to clarify the very nature of the onset at m v /2 - we 
have searched for hints of chiral symmetry restoration, both in the order parameter and 
in the spectrum, and found them to be very weak and inconclusive. The new simulations 
on large lattices have not greatly improved our understanding. We have also performed 
new, extensive simulations at strong coupling. We used different values of the bare quark 
mass, and a wide array of chemical potentials. We have also introduced new observables 
which helped shed more light on the pathologies. The analytic results available in that case 
offered a simple interpretation of the pathologies at fi = m^/2: the forbidden region of the 
simulations are to be associated with the metastable region in the effective potential. This 
suggests that the problems with finite density QCD could be solved if we could handle the 
metastable region. 

Summarizing, the simulations of quenched QCD in the scaling region, large lattices, and 
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at large \x value produce clear evidence of the major pathologies of the theory, but do not 
help to clarify the reasons behind these failures. The simulations in the strong coupling limit 
were particularly informative, especially when their results were analyzed in the context of 
available analytic treatments. 

We wish to thank I. Barbour, G. Boyd, F. Karsch, B. Petersson, E. Mendel, G. Parisi 
and F. Niedermayer for conversations and correspondence. 
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